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■ Abstract 

We present a three-dimensional numerical simulation that resolves the formation 
process of a Population III star down to a scale of ~ 100 AU. The simulation 
is initialized on the scale of a dark matter halo of mass ~ 10 6 Mq that virializes 
at z ~ 20. It then follows the formation of a fully-molecular central core, and 
traces the accretion from the diffuse dust-free cloud onto the protostellar core for 
as long as ~ 10 4 yr, at which time the core has grown to ~ 50M Q . We find that 
| the accretion rate starts very high, ~ OAMqyv' 1 , and declines rapidly thereafter 

approaching a power-law temporal scaling. Asymptotically, at times t > 10 3 yr after 
core formation, the stellar mass grows approximately as M* ~ 20MQ(t/10 3 yr) ' 4 . 
\ Earlier on, accretion is faster with M* oc t ' 75 . By extrapolating this growth over the 

full lifetime of very massive stars, t ~ 3 x 10 6 yr, we obtain the conservative upper 
limit M* < 5OOM0. The actual stellar mass is, however, likely to be significantly 
smaller than this mass limit due to radiative and mechanical feedback from the 
protostar. 
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1 Introduction 



The first (so-called Population III) stars could have had a dramatic influence 
on the early universe at redshifts z < 20. Their high yield of ionizing photons 
may have reionized the intergalactic medium (IGM) (e.g., Cen, 2003; Ciardi 
et al., 2003; Haiman and Holder, 2003; Sokasian et al., 2003; Somerville and 
Livio, 2003; Wyithe and Loeb, 2003). In addition, the first supernovae (SNe) 
were responsible for the initial IGM metal-enrichment (e.g., Madau et al., 
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2001; Mori et al., 2002; Schneider et al., 2002; Bromm et al, 2003; Furlan- 
etto and Loeb, 2003; Mackey et al., 2003; Scannapieco et al., 2003; Wada and 
Venkatesan, 2003; Norman et al., 2004; Yoshida et al., 2004). The radiative and 
hydrodynamic feedback of the first stars affected subsequent galaxy formation 
(Barkana and Loeb, 2001) and imprinted large-scale polarization anisotropies 
on the cosmic microwave background (Kaplinghat et al., 2003; Kogut et al., 
2003). In the context of popular cold dark matter (CDM) models of hier- 
archical structure formation, the first stars are predicted to have formed in 
dark matter halos of mass ~ 1O 6 M that collapsed at redshifts z ~ 20 — 30 
(e.g., Tegmark et al., 1997; Barkana and Loeb, 2001; Yoshida et al., 2003). 
Since their ionization efficiency (Tumlinson and Shull, 2000; Bromm et al., 
2001; Schaerer, 2002) and metal yield (Heger and Woosley, 2002; Umeda and 
Nomoto, 2002, 2003) depends on their mass, the fundamental question is (e.g., 
Bromm, 2004): How massive were the first stars? 

Results from recent numerical simulations of the collapse and fragmentation of 
primordial clouds suggest that the first stars were predominantly very massive, 
with typical masses M* > 100M Q (Bromm et al., 1999, 2002; Nakamura and 
Umemura, 2001; Abel et al., 2002). More specifically, these simulations have 
identified pre-stellar clumps of masses up to ~ 10 3 M Q and sizes < 0.5 pc. 
Each of the simulated clumps is conjectured to be the immediate progenitor 
of a single star or a small cluster of stars (see Larson 2003 for a review). To 
determine the mass of a single star, one must follow the fate of such a clump. 
Extending the analogous calculation for the collapse of a present-day protostar 
(Larson, 1969) to the primordial case, Omukai and Nishi (1998) have carried 
out one-dimensional hydrodynamical simulations in spherical symmetry. They 
have found that the mass of the initial hydrostatic core, formed near the center 
of the collapsing cloud when the density is sufficiently high (n ~ 10 22 cm -3 ) 
for the gas to become optically thick to continuum radiation, is almost the 
same as in present-day star formation: M core ~ 5 x 1O~ 3 M . The small value of 
the initial core has no bearing on the final mass, which is in turn determined 
by how efficient the accretion process is at incorporating the clump mass into 
the growing protostar. 

In this paper, we present results from the first three-dimensional simulation 
that is initialized on cosmological scales, followed to the formation of a proto- 
stellar high-density core, and for which the accretion flow is traced onto the 
central protostar for as long as ~ 10 4 yr after core formation. We do not yet 
self-consistently take into account the radiative and mechanical feedback from 
the protostar on the accretion flow (Omukai and Palla, 2001, 2003; Tan and 
McKee, 2003). By comparing to one-dimensional calculations that address this 
feedback on the accretion flow, we demonstrate that the accretion process is 
likely to lead to the formation of large stellar masses, as conjectured before. 

Throughout this paper, we assume a standard ACDM cosmology, with a total 
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density parameter in matter of Vt m = 1 — Vt\ = 0.3, and in baryons of Qb = 
0.045, a Hubble constant h = H /(100 km s _1 Mpc^ 1 ) = 0.7, and a power- 
spectrum amplitude <r 8 — 0.9 on 8/i _1 Mpc spheres. These values reflect the 
latest estimates of cosmological parameters from the Wilkinson Microwave 
Anisotropy Probe (WMAP) (Spergel et al, 2003). 



2 Physical ingredients 



The main physical processes operating prior to the formation of a pre-stellar 
clump with typical mass of a few hundred solar masses are described in earlier 
work (Bromm et al., 2002; hereafter BCL). Below we discuss the physical ef- 
fects that become important during the subsequent collapse to higher densities 
and the resultant formation of a protostellar core (see also Bromm, 2000). 

The formation of hydrogen molecules due to the channel saturates at a 
fractional abundance of ~ 10 -3 (BCL). At densities exceeding ~ 10 s cm~ 3 , 
however, three-body reactions are able to convert the gas into an almost fully 
molecular form (Palla et al., 1983). In addition to the reactions discussed by 
BCL (see their Table 1), we also consider here 

(13) H + H + H^H 2 + H 

(14) H + H + H 2 -> H 2 + H 2 



The reaction rates are (Palla et al., 1983): k 13 = 5.5 x 10~ 29 {T /Ky 1 cm 6 s _1 , 
and fci4 = !&i3. To estimate the density at which three-body reactions become 
important, we consider reaction (13) 

^ = M»m)' . (i) 



Defining the fractional abundance of H 2 molecules as / = 2n[u 2 ]/n-n, and 
expressing the density of hydrogen atoms n[u] = (l-f)nn with / ~ 10~ 3 1, 
we get the timescale for reaction (13) 



The value of t3_bod y equals the free-fall timescale at a critical density 



2n~~ \ 1 / 3 



f 2 Gm H 

'13 
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Evaluating £43 at T ~ 10 3 K, and assuming / ~ 10~ 3 , we find that three-body 
reactions set in at hh — 2 x 10 s cm~ 3 . Along similar lines, by considering 
reaction (14) and assuming / ~ 0.5, we find that the density at which the 
conversion of hydrogen atoms into molecules is complete is nn — 5 x 10 11 
cm -3 . Both estimates are in good agreement with our numerical results. The 
transformation of the hydrogen gas from an atomic phase to a molecular phase 
requires three modifications of the physics incorporated in the numerical code 
as described below (see also Omukai and Nishi, 1998). 

First, the heating associated with the formation of H 2 and the corresponding 
release of the molecular binding energy, en 2 = 4.48 eV, becomes important at 
high densities. To account for this additional heating source, we add to the 
energy equation (see BCL) the term (in ergs s" 1 cm" 3 ) 

r - e dn[Hal (a\ 

1 3-body — eH 2 ^ , (4 J 

for densities rin > 10 8 cm -3 . 

Second, the H2 cooling function has to be augmented to include both the 
collisional excitation of H 2 by H atoms and by H 2 molecules. The total cooling 
rate (in ergs s _1 cm -3 ) can then be written as 



A = n|/ 



IMH1 L H + /L H 2 

n H 



(5) 



where L H and L Ri (in erg s _1 cm 3 ) are the vibrational/rotational cooling co- 
efficients for collisions with H atoms and H 2 molecules, respectively. At the 
high densities considered here, all levels are close to being populated accord- 
ing to local thermodynamic equilibrium (LTE). For L H and L H2 , we use the 
parameterization given by Hollenbach and McKee (1979). 

Third, the presence of molecules leads to a modified equation of state. We 
write the gas pressure as 

P = ^—p = (7ad - l)pu , (6) 

where 7 a d is the adiabatic exponent, and u is the specific internal energy. The 
mean molecular weight, //, is given by 



where X = 0.76 and Y = 1 — X are the mass fractions in hydrogen and 
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helium, respectively. For a gas consisting of helium and pure atomic hydrogen, 
/j, ~ 1.22, whereas for a mixture of helium and pure molecular hydrogen, 
/i ~ 2.27. The adiabatic exponent can be expressed as 



(8) 




where the summation is over the 6 species H, H + , H , H 2 , e~, and He, with 
fiu 2 — 2, AiHe = 4, and /i, = 1 otherwise. In calculating the adiabatic exponent 
for H 2 , one has to consider translational, rotational, and vibrational degrees 
of freedom, resulting in 

1 ^3 + 2 + ?%) . (9) 



7h 2 - 1 2 V k B T 



Here, e v ;b is the average vibrational energy per H 2 molecule, and is given by 
(e.g., Shu, 1991) 

e 'vib = ^JZt n + , (10) 



with ftw /kB — 6300 K. For the physical conditions considered here, the 
excitation of vibrational degrees of freedom is never important. Other non- 
molecular species have only translational degrees of freedom, and consequently, 
1/(7,-1) = 3/2. We terminate the simulation when the central density is suf- 
ficiently high for opacity effects to be important. In the following, we estimate 
analytically the density beyond which the H 2 line radiation cannot escape 
from the central fully-molecular core of the gas. 

Molecular opacity starts to affect the cooling rate of the gas at a sufficiently 
high density. At temperatures T ~ 10 3 K, only the rotational transitions 
within the lowest-lying vibrational level contribute significantly to the cooling. 
Taking into account the quadrupole nature of the transitions, corresponding 
to the selection rule AJ = 2 (where J is the angular momentum quantum 
number), one may express the H 2 cooling function as 

A= Yl njAjj_ 2 AEj^ 2 . (11) 



The number density of molecules in rotational level J is given by 

nj = ^(2J+l)e-^ T , (12) 
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where n[H 2 ] is the total H 2 number density, Ej — 85.3 K J(J+ 1) is the energy- 
corresponding to level J, and Z is the partition function 



Z = ^(2J + l)e^/ feBT 



(13) 



j=o 



The energy difference between levels J and J — 2 is AEj : j_ 2 , and the cor- 
responding probability per unit time for spontaneous emission is Aj y j_ 2 oc 
AEjj_ 2 . We first identify the transition that dominates the radiative cooling 
at a given temperature, and then determine the opacity in this most important 
line. The most important transition is near the extremum set by the condition 



— {njAj y j. 2 AE Jy j_ 2 ) = , 



(14) 



implying 



j ~ fi 4 [ 

max VlOOOK 



(15) 



where T ~ 10 3 K is typical for the central region of the collapsing clump. 



For the 6^4 transition, the line absorption coefficient (in cm 1 ) can be 
expressed as (Rybicki and Light man, 1979) 



AE, 



OL v 



6,1 



47T 



n 4 B 4fi 



I _ e -AE 6A /k B T 



(16) 



where we have taken the levels to be populated according to LTE. By assum- 
ing that broadening is entirely due to the thermal Doppler effect, the profile 
function is 



y/7rAu D 



cxp 



Avl 



(17) 



with Avd = (z/ /c)y 2/cBT/m, a particle mass m = Im^ for H 2 , and /iz/ = 
AE 6i4 . This function leads to a modest overestimate of the opacity, as parts 
of the accretion flow are mildly supersonic (with a maximum Mach number of 
~ 1.5 at radius ~ 200 AU), but our goal here is only to get a rough estimate 
as to when line opacity becomes important. In the following, we approximate 
tp(v) ~ <p(i/ ) = l/(v^Az/ D ), and use Au D /u ~ 10~ 5 for T ~ 10 3 K. The 
Einstein coefficient for stimulated absorption is 

B « = ? 2?ot- 4w * °- 04 cm2 erg " s " ■ (18) 
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where A§^ is taken from Turner et al. (1977). At T ~ 10 3 K, the partition 
function is Z ~ 12.1, and using equation (15) we find n 4 ~ 0.1n[H 2 ]. 

We may now estimate the density at which the optical depth, r v ~ a u L, 
approaches unity, where L is the characteristic size of the high-density, fully 
molecular region 



From our simulation, we estimate L ~ 10~ 4 pc, and consequently have riu ~ 
10 11 cm -3 . At densities exceeding this value, radiation can still escape in less 
opaque lines with a smaller A coefficient, and in the continuum between the 
lines, but neglecting the effects of radiative transfer renders the results of the 
simulation increasingly unreliable at yet higher densities. Omukai and Nishi 
(1998), who did include the treatment of radiative transfer, were able to follow 
the collapse all the way to stellar densities, although their approach was limited 
to a spherical geometry. In our numerical simulation, we form a sink particle 
once the gas density exceeds a value of 10 12 cm -3 (see BCL for details about 
the numerical implementation). This sink particle approximates a hydrostatic 
core at the center of the accretion flow, although it has a size that is still 
much larger than the expected scale of the actual core (Omukai and Nishi, 
1998; Ripamonti et al., 2002). 



3 Numerical methodology 

The evolution of the dark matter and gas components is calculated with a 
version of TREESPH (Hernquist and Katz, 1989), combining the smoothed 
particle hydrodynamics (SPH) method with a hierarchical (tree) gravity solver 
(see BCL for further details). Here, we briefly describe the modifications in- 
troduced to the code in order to allow it to treat accretion flows in a dense, 
high-redshift environment. 

We initialize the simulation (see § 4.1) on the scale of the minihalo (> 150 pc), 
and follow the collapse to the point where the gas is converted into a fully 
molecular phase (on a scale < 100 AU). In order to treat the large dynamic 
range spanned between the cosmological and protostellar length scales, we 
carry out our simulation in two steps. We first simulate the collapse of the 
primordial gas, together with the virialization of the dark matter (DM) in 
the minihalo, until the formation of a gravitationally-bound clump (similar to 
BCL). At this point, we stop the simulation, and resample the central density 
field (involving ~ 3000M Q in gas within a spherical region of radius 3.1 pc), 




(19) 
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with an increased number of SPH particles. We adopt a rapidly decreasing 
timestep according to the scaling At sys oc l/^/Gp max . Here, p max is the max- 
imum gas density at a given instant and the system timestep, At sys , is the 
maximum time by which a particle is advanced within the multiple-timestep 
scheme employed by the simulations (see Hernquist and Katz, 1989, for de- 
tails). Over the short timestep dictated by the maximum density, the global 
system does not evolve by much, while the runaway collapse of the compact 
dense core converges rapidly. 

In general, the mass resolution of an SPH simulation is approximately 



where iV ne igh — 32 is the number of neighboring particles within a given SPH 
smoothing kernel, iVspH the total number of SPH particles, and Mb the total 
baryonic mass. To avoid numerical artifacts the Jeans mass, Mj, has to be 
resolved (Bate et al., 2003), and so we require M res < Mj. After the onset 
of gravitational instability in the large-scale simulation, the rapid increase in 
density naturally would lead to a violation of this criterion. The resampling 
technique, however, allows us to follow the collapse to a much higher density 
due to the improved M res . 

We here briefly summarize the resampling procedure (see Bromm, 2000, and 
Bromm and Loeb, 2003, for details). Every SPH particle in the original, un- 
refined simulation acts as a parent particle p, and spawns iV re f child particles 
i, where iV re f = 50 in this paper. The child particles are distributed according 
to the SPH smoothing kernel W{fi — r p ; h p ), by employing a standard Monte- 
Carlo comparison-rejection method. Here, h p is the smoothing length of the 
parent particle. The velocity, temperature, H2 abundance, and free-electron 
abundance are directly inherited from the parent particle, = v p , Tj = T p , 
fi — fpi an d Xi = x p for all i, respectively. Finally, each child particle is as- 
signed a mass rrii = m p /N ref . This procedure guarantees the conservation of 
linear and angular momentum. Energy is also conserved well, although there is 
a small artificial contribution to the gravitational potential energy due to the 
discreteness of the resampling. The resampling results in iVgpH = 35,500 within 
the central high- resolution region, and the mass resolution is now M res — 4M , 
as compared to M rcs ~ 2OOM in the original simulation. Our technique of 
refining a coarser, parent simulation, and following the further collapse with 
increased resolution, is conceptually similar to the adaptive mesh refinement 
(AMR) method which has already been successfully applied to problems in 
star formation (e.g., Truelove et al., 1998; Abel et al., 2002), but our approach 
is one of the first attempts of implementing such a scheme within SPH (see 
also Kitsionas and Whitworth, 2002; Bromm and Loeb, 2003). 




(20) 
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To study accretion onto the central core, the diffuse gas in the envelope must 
be followed for a sufficiently long time (roughly equal to the local free-fall 
timescale). It is, therefore, necessary to introduce a central sink particle at 
some stage in the evolution. Otherwise, we would have to adopt smaller and 
smaller timesteps to be able to follow the central collapse, and with limited 
computer resources the surrounding envelope would not have time to be ac- 
creted onto the core. BCL prescribed SPH particles to merge into more massive 
ones, provided they exceed a pre-determined density threshold. The sink par- 
ticle technique allows us to study the ongoing accretion process over many 
dynamical times. 



4 Simulating the accretion flow 

4-1 Initial Conditions 

Within the hierarchical ACDM model, the first luminous objects are expected 
to form out of high-cx peaks in the random field of primordial density fluctua- 
tions. The early (linear) evolution of such a peak, assumed to be isolated and 
roughly spherical, can be described by the top-hat model (e.g., Padmanab- 
han, 1993). We use the top-hat approximation to specify the initial amplitude 
of the CDM configuration. In this paper, we investigate the fate of a high-a 
peak of total mass 10 6 M Q , corresponding to 1.5 x 1O 5 M in baryons, which 
collapses (or virializes) at a redshift z vil ~ 20. 

Our simulation is initialized at Zi = 100, by performing the following steps. 
The collisionless DM particles are placed on a cubical Cartesian grid, and are 
then perturbed according to a given power spectrum P(k) = Ak n , by apply- 
ing the Zeldovich (1970) approximation which also allows to self-consistently 
assign initial velocities. The power-law index is set to n = — 3 which approxi- 
mately describes the spectral behavior on the mass scale of ~ 1O 6 M . To fix 
the amplitude A, we specify the initial variance of the fluctuations 

aj = A^k n . (21) 

The summation is over all contributing modes, where the minimum wavenum- 
ber is given by the overall size of the Cartesian box, and the maximum 
wavenumber is set by the Nyquist frequency, k max ~ 2n/d gr id- Here, <i gri d is the 
mesh size of the Cartesian grid. Choosing Oi ~ 0.2 at Zi = 100, the rms fluc- 
tuation at the moment of collapse is a(z = 20) = [(1 + Zj) /(l + z)]o~i ~ 1. This 
choice ensures that substructure develops on a similar timescale as the overall 
collapse of the background medium. Next, particles within a (proper) radius of 
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Fig. 1. Collapse and fragmentation of a primordial cloud. We show the projected 
gas density at a redshift z — 21.5, shortly after gravitational runaway collapse has 
converged at the center of the cloud. Left: The coarse-grained morphology in a box 
with a linear physical size of 23.5 pc. At this time, a gravitationally-bound clump has 
formed with a mass of ~ 10 3 Mq. Right: The fine-grained morphology in a box with 
a linear physical size of 0.5 pc. The central density peak accretes mass vigorously, 
and is accompanied by a secondary clump. 

Ri ~ 160 pc are selected for the simulation. The resulting number of DM par- 
ticles is Ndm = 17, 074. Finally, the particles are set into rigid rotation and are 
endowed with a uniform Hubble expansion (see Katz, 1991). Angular momen- 
tum is added by assuming a spin parameter A s = L\E\ 1/2 / (GM 5/2 ) = 0.05, 
close to the mean- value measured in cosmological N-body simulations (e.g., 
Jang-Condell and Hernquist, 2001). Here, L, E, and M are the total angular 
momentum, binding energy, and mass, respectively. The spin parameter is a 
measure of the degree of rotational support, such that the ratio of centrifugal 
to gravitational acceleration is given by ~ at virialization. The collisional 
SPH particles (Asp H = 32, 768) are randomly placed to approximate a uni- 
form initial density. The SPH particles follow the same Hubble expansion and 
rigid rotation as the DM particles. 

4-2 Results 

The initial evolution up to the formation of a pre-steliar clump is very similar 
to the cases studied by BCL. The gas dissipatively settles into the center of 
the CDM minihalo, reaching the characteristic state with density ~ 10 4 cm -3 
and temperature ~ 200 K. In Figure 1 (left panel), we show the morphology 
within the central ~ 25 pc, briefly after the first high-density core has formed 
clS cL result of gravitational runaway collapse. 

To study the formation of a Population III star, we have re-simulated the evo- 
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Fig. 2. Radial distribution of physical properties of SPH particles within the main 
clump. The presented snapshot was taken shortly before the central sink particle was 
created, (a) Hydrogen number density (in cm~ 3 ) vs. radial distance from density 
maximum (in AU). For comparison, we also show the isothermal profile (dashed 
line), which provides a good description for the outer region, (b) Gas temperature 
(in K). (c) Radial velocity (in km s _1 ). Dashed line: Radially averaged (negative) 
sound speed. It is evident that the flow is mostly subsonic, with a few fluid elements 
having mildly supersonic speeds, (d) Hydrogen molecule fraction. Within the central 
~ 100 AU, three-body reactions have converted the gas into a fully molecular phase. 



lution of the central clump with sufficient resolution to follow its collapse up 
to a limiting density of n ~ 10 12 cm -3 , at which point opacity effects become 
important (see §2), and a sink particle is created. The right panel in Figure 1 
shows the gas density on a scale of 0.5 pc. Several features are evident in this 
image. First, the central clump does not undergo further sub-fragmentation, 
and is likely to form a single Population III star. Second, a companion clump 
is visible at a distance of ~ 0.25 pc. If negative feedback from the first-forming 
star is ignored, this companion clump would undergo runaway collapse on its 
own, approximately ~ 3 Myr later. This timescale is comparable to the life- 
time of a very massive star (e.g., Bromm et al., 2001). If the second clump 
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were able to survive the intense radiative heating from its neighbor, it could 
become a star before the first one explodes as a SN. Whether more than one 
star can form in a low-mass halo thus crucially depends on the degree of syn- 
chronization of clump formation. Finally, the non-axisymmetric disturbance 
induced by the companion clump, as well as the angular momentum stored 
in the orbital motion of the binary system, allow the system to overcome the 
angular momentum barrier for the collapse of the central clump (see Larson, 
2002), and to possibly form a wide binary star system. 

In Figure 2, we consider the radial structure around the central density maxi- 
mum, shortly before the accreting sink particle is created. The overall behavior 
is similar to previous results (Omukai and Nishi, 1998; Abel et al., 2002; Ri- 
pamonti et al., 2002). The density profile [panel (a)] consists of a central, flat 
core, surrounded by a roughly isothermal envelope. Such a configuration is 
generically predicted for collapsing star forming clumps (see Larson, 2003), 
and can approximately be described by the Larson-Penston similarity solu- 
tion. Three-body reactions have succeeded in converting the central > 1M Q 
into fully molecular form [panel (d)]. Notice that the central temperature never 
drops precipitously due to the boost in H 2 cooling, even though H 2 is now more 
abundant by a factor of ~ 10 3 . The opposite effects of compressional heating 
and H 2 formation heating balance the enhanced cooling rate, and lead to an 
almost isothermal collapse in the central region. 

Next, we address the accretion flow onto the central sink particle by tracing 
its growth for ~ 10 4 yr after its initial formation. This accretion process will 
determine how massive the star finally gets. We begin with a brief discussion 
of the basic physics of the accretion, then describe the results from earlier ID 
calculations, and finally compare these results to our simulation. In general, 
star formation typically proceeds from the 'inside-out' through the accretion 
of gas onto a central hydrostatic core (e.g., Larson, 2003). Even though the 
initial mass of the hydrostatic core in the Population III case is very similar 
to that in the Population I/II case, the accretion process is expected to be 
rather different. On dimensional grounds, the accretion rate simply scales as 
the sound speed cubed over Newton's constant (or equivalently the ratio of 
the Jeans mass and the free-fall time): M acc ~ c 3 s /G oc T 3 / 2 . A comparison of 
the temperatures in present-day star forming regions (T ~ 10 K) with those 
in primordial ones (T ~ 200 — 300 K) already indicates a difference in the 
accretion rate of more than two orders of magnitude. 

Omukai and Palla (2001, 2003) extended earlier work by Stahler et al. (1986) 
and investigated the spherical accretion problem in considerable detail, go- 
ing beyond the simple dimensional argument given above. Their computa- 
tional technique approximates the time evolution by considering a sequence 
of steady-state accretion flows onto a growing hydrostatic core. Somewhat 
counterintuitively, these authors identified a critical accretion rate, M crit ~ 
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4 x 10~ 3 M Q yr -1 , beyond which the protostar cannot grow to masses much in 
excess of a few tens of solar masses. The critical rate is predicted to be almost 
independent of mass (Omukai and Palla, 2003). For smaller rates, however, 
the accretion is predicted to proceed all the way up to ~ 6OOM , i.e., of or- 
der the host clump mass. The physical origin of the critical accretion rate is 
that for ongoing accretion onto the core, the luminosity must not exceed the 
Eddington value, Ledd, if the inner region of the gaseous core is ionized by 
the central star. In our simulation, we do not resolve the evolution of the fully 
molecular core to temperatures and densities high enough to completely ion- 
ize the gas. Our estimate of the accretion rate, evaluated at the density where 
the gas becomes fully molecular, is likely to be applicable also at the higher 
density where the gas becomes fully ionized. This extrapolation assumes that 
the inflowing gas will not be diverted into outflows before reaching the cen- 
tral protostar. It will be important to test this assumption in future work, 
taking into account the magneto-hydrodynamics of possible bipolar jets. We 
now ask whether the growing star is able to settle onto the main sequence 
(MS), without experiencing a strong radiative feedback on the accretion flow. 
Before the onset of hydrogen burning, the luminosity is approximately given 
by L tot ~ L acc ~ GM*M acc /-R*. Here, we have ignored the contribution to the 
overall luminosity from Kelvin-Helmholtz (KH) contraction for simplicity, as 
we are only aiming at an order-of-magnitude estimate. By setting L acc ~ Ledd 
it follows that 

M crit ^ ~ 5 x 1O- 3 M yr- 1 , (22) 



where /?* ~ 5R & , a typical value for a Population III MS star (e.g., Bromm et 
al., 2001). This is the critical accretion rate that must not be exceeded upon 
reaching the MS. If it is surpassed, radiation pressure on free electrons will 
drive a strong wind, thus preventing further accretion. 

It is possible to start out with accretion rates that are significantly larger 
than M cri t, because at earlier times the protostellar radius is much larger, 
R* > 100-R©, and it only gradually shrinks to the MS value in the course 
of KH contraction. During the KH phase, the star is in hydrostatic equilib- 
rium, but not yet in thermal equilibrium, which is only achieved on the MS 
when contraction temporarily stops. The above condition can only determine 
whether accretion is capable of incorporating most of the diffuse envelope into 
the star, or whether it is shut-off by radiation pressure early on. The precise 
stellar mass in the latter case can only be determined by a detailed simulation 
that takes both radiative-transfer and the hydrodynamics into account. Such a 
fully self-consistent calculation has not yet been achieved in three dimensions, 
although Omukai and Palla (2003), as well as Tan and McKee (2003), have 
treated this demanding problem in zero- and one dimension. 
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Fig. 3. Accretion onto a primordial protostar. The morphology of the accretion 
flow is shown in Fig. 1. Left: Accretion rate (in M Q yr -1 ) vs. time (in yr) after 
the formation of a molecular core. Solid line: Accretion rates approximated as: 
M acc oc i _a25 for t < 10 3 yr, and M acc oc t~ oe for t > 10 3 yr. The early 3D evolution 
does not follow the simple power-law temporal behavior found in ID simulations 
with very high resolution (Omukai and Nishi, 1998). Right: Mass of the central core 
(in Mq) vs. time. Solid lines: Accretion history approximated as: oc t 0,75 for 
t < 10 3 yr, and M* oc t 0A for t > 10 3 yr. Using this analytical approximation, wc 
extrapolate that the protostellar mass grows to ~ 12OM after ~ 10 5 yr, and to 
~ 5OOM after ~ 3 x 10 6 yr, the lifetime of a very massive star. 



Realistic flows are expected to produce a time-dependent accretion rate, and so 
their outcome depends on whether the accretion rate declines rapidly enough 
to avoid exceeding the Eddington luminosity at some stage during the evolu- 
tion. The main caveat concerning the Omukai and Palla (2001, 2003) results 
involves the geometry of the infall. A three-dimensional accretion flow of gas 
with angular momentum will deviate from spherical symmetry, and instead 
form a flattened rotating configuration such as a disk. In this case, most of the 
photons will likely escape along the rotation axis (where the gas density and 
the corresponding opacity is low) whereas mass may flow unimpeded along 
the disk plane (see Tan and McKee, 2003). 

In Figure 3, we show the accretion history. Our time-dependent, 3-dimensional 
simulation provides the accretion flow around a primordial protostar and 
shows how the molecular core grows in mass over the first ~ 10 4 yr after 
its formation. The accretion rate is initially very high, M acc ~ O.1M yr^ 1 , 
and subsequently declines with time in a complex fashion. Expressed in terms 
of the sound speed in Population III pre-stellar cores, the initial rate is: 
M acc ~ 25c 3 /G. Initial accretion rates of a few tens times c 3 /G are commonly 
encountered in simulations of present-day star formation (Larson, 2003). We 
thus find that Population III stars reproduce this generic behavior. The mass 
of the molecular core, taken as an estimate for the protostellar mass, grows 
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approximately as 

/ t \ °' 75 

M* = J M acc dt » 2OM I j (23) 



for t < 10 3 yr, and as 

M* = J M acc dt « 2OM ( ^— J (24) 



afterwards. The early-time accretion behavior is similar to the power-law scal- 
ings found in ID, spherically symmetric simulations (Omukai and Nishi, 1998; 
Ripamonti et al., 2002). The complex accretion history encountered in our 
simulation, however, cannot be accommodated with a single power law. Such 
a multiple power-law behavior has previously been inferred (Abel et al., 2002; 
Omukai and Palla, 2003). This earlier inference, though, was indirect, in that 
it did not follow the actual accretion flow, but instead estimated the accretion 
from the instantaneous density and velocity profiles at the termination of the 
3D simulation of Abel et al. (2002). Expressed as a function of stellar mass, 
our accretion rate scales as M acc oc M~ 0-35 for t < 10 3 yr, and M acc oc M~ 1-5 
afterwards. The early behavior is similar to the accretion law predicted by 
Tan and McKee (2003): M acc oc M~ 0A . In summary, our simulated accretion 
history at early times is consistent with previous estimates, based on 0D and 
ID work, but deviates significantly at later times, at t > 10 3 yr. Further- 
more, the accretion flow is not a simple scale-free process when followed over 
a sufficiently long time. The self-similarity of the flow is broken because the 
presence of three-body reactions introduces a preferred timescale (see equ. 
(2)): t 3 -body ~ 10 3 yr. 

Our simulations do not take into account any feedback effects from the proto- 
star on the accretion flow. To ascertain when these are likely to intervene, we 
compare our accretion history with the critical threshold derived by Omukai 
and Palla (2003). We find that as early as ~ 3000 yr after the accretion 
started, the accretion rate drops below M cr it, at which point the stellar mass 
is M* ~ 3OM . The accretion rate becomes sub-critical early on, while the 
protostar still undergoes KH contraction. Thus, according to the model of 
Omukai and Palla (2003), the star can settle onto the MS without experienc- 
ing a strong continuum-driven wind, and accretion can continue during the 
hydrogen-burning stage. A conservative upper limit for the final mass of the 
star is then, M*(t = 3 x 10 6 yr) ~ 5OOM , assuming that accretion cannot go 
on for longer than the total lifetime of a very massive star, which is almost 
independent of stellar mass (e.g., Bond et al., 1984; Bromm et al., 2001). 



15 



5 Conclusions 



We have presented the first three-dimensional simulation of Population III 
star formation that followed the accretion flow for ~ 10 4 yr, allowing to better 
estimate the final stellar mass. The earlier collapse phase prior to the forma- 
tion of the initial hydrostatic core, gives results that are similar to previous 
work (Omukai and Nishi, 1998; Abel et al, 2002; Ripamonti et al., 2002). In 
particular, we confirm that the gas is converted into a fully molecular phase 
due to three-body reactions within the central ~ 100 AU (comprising an ini- 
tial mass of > 1M ). When extrapolating the accretion to ~ 3 x 10 6 yr, the 
lifetime of a massive star, we get M* ~ 5OOM . 

Can a Population III star ever reach this asymptotic mass limit? Our hydro- 
dynamic simulation lacks radiative transfer and cannot address this question, 
as the answer depends on whether the accretion from the dust-free envelope 
is eventually terminated by feedback from the star (e.g., Omukai and Palla, 
2001, 2003; Ripamonti et al., 2002; Omukai and Inutsuka, 2002; Tan and 
McKee, 2003). The standard mechanism by which accretion may be termi- 
nated in metal-rich gas, namely radiation pressure on dust grains (Wolfire 
and Cassinelli, 1987), is obviously not effective for gas with a primordial com- 
position. It has been speculated that accretion could instead be turned off 
through the formation of an H II region (Omukai and Inutsuka, 2002), or 
through the radiation pressure exerted by trapped Lya photons (Tan and 
McKee, 2003). The unsolved problem of whether the accretion process is ter- 
minated by radiative feedback from the protostar, defines the current frontier 
in three-dimensional numerical studies of the formation of Population III stars. 
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